Load libraries
source("../General/parameters.R")
source("../General/rmb_functions.R")
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2.9000 ✓ purrr 0.3.4
✓ tibble 3.0.1 ✓ dplyr 1.0.0
✓ tidyr 1.1.0 ✓ stringr 1.4.0
✓ readr 1.3.1 ✓ forcats 0.5.0
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(cowplot)
********************************************************
Note: As of version 1.0.0, cowplot does not change the
default ggplot2 theme anymore. To recover the previous
behavior, execute:
theme_set(theme_cowplot())
********************************************************
library(randomForest)
randomForest 4.6-14
Type rfNews() to see new features/changes/bug fixes.
Attaching package: ‘randomForest’
The following object is masked from ‘package:dplyr’:
combine
The following object is masked from ‘package:ggplot2’:
margin
library(tidyverse)
Load data
map <- readRDS("../Data/drought_map.RDS")
tax <- readRDS("../Data/tax.RDS")
otu <- readRDS("../Data/otu_pers.RDS")
train.set <- readRDS("../Data/rf_results.RDS") %>%
mutate(Compartment = fct_recode(Compartment,
"Rhizosphere" = "RS",
"Endosphere" = "ES")) %>%
mutate(Compartment = fct_relevel(Compartment, "Rhizosphere"))
all.clust.summary <- readRDS("../Data/drought_clusters.RDS")
Load data
## Get log-transformed relative abundances for the OTU table
otu <- rel_ab(otu)
otu <- log_norm(otu)
map <- mutate(map, Set = ifelse(Treatment2 == "WC_TRN", "Train", "Test"))
# The predictors need to be formatted as a data frame or matrix
get_predictors <- function(x){
select(x, -(SampleID:Age))
}
# The response variable needs to be formatted as a vector
get_response <- function(x) {
select(x, Age) %>%
.$Age
}
whole.set <- as.data.frame(t(otu)) %>%
mutate(SampleID = rownames(.)) %>%
inner_join(select(map, SampleID, Compartment, Set, Treatment, Age), by = "SampleID") %>%
group_by(Set, Compartment, Treatment) %>%
nest() %>%
mutate(otu = map(data, get_predictors),
age = map(data, get_response)) %>%
mutate(Compartment = fct_recode(Compartment,
"Rhizosphere" = "RS",
"Endosphere" = "ES")) %>%
mutate(Compartment = fct_relevel(Compartment, "Rhizosphere"))
trt.lines <- data.frame(Treatment = c("DS1", "DS2", "DS3"),
Treatment2 = c("DS1", "DS2", "DS3"),
Contrast = c("WC vs DS1", "WC vs DS2", "WC vs DS3"),
IniTreatment = c(4.5, 4.5, 4.5),
EndTreatment = c(5.5, 6.5, 7.5),
IniTreatment2 = c(41,41,41),
EndTreatment2 = c(52,62,74))
cv.plot <- train.set %>%
unnest(cv, nOTU) %>%
mutate(Compartment = fct_relevel(Compartment, "RS")) %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
ggplot(aes(nOTU, Error, color = Compartment)) +
geom_line(size =1) +
geom_point(aes(alpha = minError), size = 10) +
geom_text(aes(label = nOTU, alpha = minError), size = 6, color = "black") +
scale_color_manual(values = rev(cmp.pal)) +
scale_alpha_manual(values = c(0,1)) +
guides(alpha = F) +
scale_x_log10() +
theme_classic() +
theme(text = element_text(size = 15),
legend.position = c(0.5, 0.8),
legend.title = element_blank())
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(cv, nOTU))`, with `mutate()` if neededUnknown levels in `f`: RSUnknown levels in `f`: RS, ES
cv.plot
Check taxonomy of the OTUs
rf.ranks <- train.set %>%
unnest(imp.tax) %>%
mutate(cmpOTU = paste(Compartment, OTU_ID),
Rank = rank(PercIncMSE))
rf.ranks %>%
ggplot(aes(reorder(cmpOTU, Rank), PercIncMSE, fill = PhyClass2)) +
geom_bar(stat = "identity") +
xlab("OTU") +
scale_y_log10() +
coord_flip() +
scale_fill_manual(values = phy.pal) +
facet_wrap(~Compartment, scales = "free") +
theme_light() +
theme(text = element_text(size = 15),
axis.text.y = element_blank())
otu.rf <- train.set %>%
unnest(imp.tax) %>%
mutate(cmpOTU = paste(Compartment, OTU_ID),
Rank = rank(PercIncMSE))
rs.id <- filter(otu.rf, Compartment == "Rhizosphere")$OTU_ID
es.id <- filter(otu.rf, Compartment == "Endosphere")$OTU_ID
# otu2 <- readRDS("../Data/otu_pers.RDS")
# otu2 <- rel_ab(otu2)
# otu2 <- log_norm(otu2)
otu.tidy <- tidy_otu(otu)
get_order <- function(df, var) {
mtrx <- df %>%
mutate(Group = paste(Compartment, Treatment, Time, sep = ".")) %>%
select(Group, var, OTU_ID) %>%
spread(key = Group, value = var)
mtrx <- as.data.frame(mtrx)
rownames(mtrx) <- mtrx$OTU_ID
mtrx <- mtrx[,-1]
dist.tmp <- dist(as.matrix(mtrx))
dist.clust <- hclust(dist.tmp, method = "average")
ord.names <- dist.clust$labels[dist.clust$order]
clust.ord <- data.frame(OTU_ID = ord.names, order = 1:length(ord.names))
df <- inner_join(df, clust.ord, by = "OTU_ID")
df$order
}
get_cluster <- function(df, var, nclust) {
mtrx <- df %>%
mutate(Group = paste(Compartment, Treatment, Time, sep = ".")) %>%
select(Group, var, OTU_ID) %>%
spread(key = Group, value = var)
mtrx <- as.data.frame(mtrx)
rownames(mtrx) <- mtrx$OTU_ID
mtrx <- mtrx[,-1]
dist.tmp <- dist(as.matrix(mtrx))
dist.clust <- hclust(dist.tmp, method = "average")
ord.names <- dist.clust$labels[dist.clust$order]
clust.ord <- data.frame(OTU_ID = ord.names, order = 1:length(ord.names))
clust.cut <- cutree(dist.clust[c(1,2,4)], k = nclust)
clust.ord$Cluster <- as.factor(clust.cut[clust.ord$OTU_ID])
df <- inner_join(df, clust.ord, by = "OTU_ID")
df$Cluster
}
es.master.tidy <- otu.tidy %>%
inner_join(select(map, SampleID, Time, Compartment, Treatment), by = "SampleID") %>%
filter(Compartment == "ES") %>%
group_by(Compartment, Treatment, OTU_ID) %>%
mutate(zscore = (Count - mean(Count))/sd(Count)) %>%
group_by(Compartment, Treatment, Time, OTU_ID) %>%
summarise(MeanZS = mean(zscore)) %>%
ungroup()
`summarise()` regrouping output by 'Compartment', 'Treatment', 'Time' (override with `.groups` argument)
es.order <- es.master.tidy %>%
filter(OTU_ID %in% es.id) %>%
filter(Treatment == "WC") %>%
mutate(ZS_order = get_order(., "MeanZS"),
ZS_cluster = get_cluster(., "MeanZS", 3)) %>%
mutate(Trend = case_when(
ZS_cluster == 1 ~ "Complex",
ZS_cluster == 2 ~ "Early Colonizer",
ZS_cluster == 3 ~ "Late Colonizer"
))
Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(var)` instead of `var` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.
es.plot.df <- es.master.tidy %>%
inner_join(select(es.order, OTU_ID, ZS_order, ZS_cluster, Trend), by = "OTU_ID") %>%
mutate(Treatment = fct_relevel(Treatment, "WC")) %>%
inner_join(tax, by = "OTU_ID") %>%
mutate(EndTreatment = case_when(
Treatment == "WC" ~ 4.5,
Treatment == "DS10" ~ 5.5,
Treatment == "DS20" ~ 6.5,
Treatment == "DS30" ~ 7.5,
))
rs.master.tidy <- otu.tidy %>%
inner_join(select(map, SampleID, Time, Compartment, Treatment), by = "SampleID") %>%
filter(Compartment == "RS") %>%
group_by(Compartment, Treatment, OTU_ID) %>%
mutate(zscore = (Count - mean(Count))/sd(Count)) %>%
group_by(Compartment, Treatment, Time, OTU_ID) %>%
summarise(MeanZS = mean(zscore)) %>%
ungroup()
`summarise()` regrouping output by 'Compartment', 'Treatment', 'Time' (override with `.groups` argument)
rs.order <- rs.master.tidy %>%
filter(OTU_ID %in% rs.id) %>%
filter(Treatment == "WC") %>%
mutate(ZS_order = get_order(., "MeanZS"),
ZS_cluster = get_cluster(., "MeanZS", 3)) %>%
mutate(Trend = case_when(
ZS_cluster == 3 ~ "Complex",
ZS_cluster == 1 ~ "Early Colonizer",
ZS_cluster == 2 ~ "Late Colonizer"
))
rs.plot.df <- rs.master.tidy %>%
inner_join(select(rs.order, OTU_ID, ZS_order, ZS_cluster, Trend), by = "OTU_ID") %>%
mutate(Treatment = fct_relevel(Treatment, "WC")) %>%
inner_join(tax, by = "OTU_ID") %>%
mutate(EndTreatment = case_when(
Treatment == "WC" ~ 4.5,
Treatment == "DS10" ~ 5.5,
Treatment == "DS20" ~ 6.5,
Treatment == "DS30" ~ 7.5,
))
all.plot.df <- rbind(rs.plot.df, es.plot.df) %>%
mutate(OTU_ID = paste(Compartment, OTU_ID))
all.ord <- rbind(group_by(rs.order, Compartment, OTU_ID, ZS_cluster, Trend) %>% count(),
group_by(es.order, Compartment, OTU_ID, ZS_cluster, Trend) %>% count())
all.ord$Trend <- as.factor(all.ord$Trend)
trend.p <- all.plot.df %>%
mutate(Treatment = fct_recode(Treatment, WC = "WC", DS1 = "D1", DS2 = "D2", DS3 = "D3")) %>%
mutate(Compartment = fct_relevel(Compartment, "RS")) %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
mutate(Trend = fct_relevel(Trend, "Complex", after = Inf)) %>%
ggplot(aes(Time * 10, reorder(OTU_ID, ZS_order), fill = MeanZS)) +
geom_tile() +
geom_point(aes(x = -1, color = Trend)) +
geom_vline(data = trt.lines, aes(xintercept = (IniTreatment * 10)), linetype = 3) +
geom_vline(data = trt.lines, aes(xintercept = (EndTreatment* 10)), linetype = 3) +
scale_fill_gradientn(colors = RColorBrewer::brewer.pal(9, "Greys"),
name = "Relative Abundance\n(z-score)") +
scale_color_manual(values = c("dodgerblue", "gold", "grey50")) +
ylab("Age Discriminant OTU") +
xlab("Plant Age (days)") +
theme_classic() +
theme(axis.text.y = element_blank(),
panel.grid = element_blank(),
text = element_text(size = 15),
legend.position = "bottom",
axis.line.y = element_blank(),
axis.ticks.y = element_blank()) +
facet_grid(Compartment ~ Treatment, scales = "free", space = "free") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5))
trend.p
rf.otus <- rbind(rs.plot.df, es.plot.df) %>%
group_by(Compartment, OTU_ID, Trend, Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
count() %>%
select(-n) %>%
ungroup()
write.table(rf.otus, "../Tables/rf_otus.tsv", sep = "\t", quote = F, row.names = F)
overlap.plot <- all.plot.df %>%
separate(OTU_ID, c("tmp", "OTU_ID")) %>%
select(OTU_ID, Compartment, Trend) %>%
mutate(RFTrend = Trend) %>%
select(-Trend) %>%
group_by(OTU_ID, Compartment, RFTrend) %>%
count() %>%
select(-n) %>%
left_join(all.clust.summary, by = c("Compartment", "OTU_ID")) %>%
ungroup() %>%
mutate(Compartment = fct_relevel(Compartment, "RS")) %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
mutate(Trend = ifelse(is.na(Trend), "None", Trend)) %>%
mutate(Trend = fct_recode(Trend,
"Semi-Persistent\nEnrichment" = "Semi-Persistent Enrichment",
"Persistent\nDepletion" = "Persistent Depletion",
"Transient\nDepletion" = "Transient Depletion",
"Transient\nEnrichment" = "Transient Enrichment")) %>%
mutate(Trend = fct_relevel(Trend, "None", after = Inf)) %>%
mutate(RFTrend = fct_relevel(RFTrend, "Early Colonizer", after = Inf)) %>%
ggplot(aes(RFTrend, fill = Trend)) +
geom_bar() +
facet_grid(Compartment ~.) +
ylab("nOTU") +
xlab("") +
scale_fill_manual(name = "Drought Module",
values = c(RColorBrewer::brewer.pal(4,"Dark2"), "gray69")) +
guides(fill = guide_legend(ncol = 1)) +
theme_classic() +
theme(text = element_text(size = 15),
legend.position = "right") +
coord_flip()
overlap.plot
rf.tax <- all.plot.df %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
mutate(Compartment = fct_relevel(Compartment, "Rhizosphere")) %>%
mutate(Trend = fct_relevel(Trend, "Complex", after = Inf)) %>%
mutate(Trend = fct_recode(Trend, Early = "Early Colonizer", Late = "Late Colonizer")) %>%
group_by(Compartment, Trend, OTU_ID) %>%
count() %>%
separate(OTU_ID, c("tmp", "OTU_ID")) %>%
select(-n, -tmp) %>%
inner_join(tax, by = "OTU_ID") %>%
group_by(Compartment, Trend, PhyClass2, Phylum, Class, Order) %>%
count() %>%
ggplot(aes(Trend, paste(Phylum, Class, Order, sep = " / "), fill = PhyClass2, label = n)) +
geom_tile(color = "white", size = 1) +
geom_text(color = "black", size = 3) +
scale_fill_manual(name = "Taxon",
values = phy.pal,
limits = levels(tax$PhyClass2)) +
ylab("Phylum / Class / Order") +
facet_grid(.~Compartment) +
theme_classic() +
theme(text = element_text(size = 15),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_text(size = 10),
legend.position = "none",
panel.grid.major = element_line(colour = "gray90"),
axis.line.y = element_blank(),
axis.ticks.y = element_blank())
rf.tax
supp.left <- plot_grid(cv.plot,
overlap.plot,
nrow = 1,
rel_widths = c(2,3),
labels = c("A", "B"),
label_size = 20)
supp.left
plot_grid(supp.left,
trend.p,
ncol = 1,
rel_heights = c(1,3),
labels = c(NA, "C"),
label_size = 20)
Run random forest using the top OTUs as predictors
get_rf <- function(imp.predictors, response) {
randomForest(imp.predictors, response, importance = F, keep.forest = T)
}
train.set <- train.set %>%
mutate(rf = map2(imp.otu, age, get_rf))
train.set %>%
filter(Compartment == "Rhizosphere") %>%
.$rf
[[1]]
Call:
randomForest(x = imp.predictors, y = response, importance = F, keep.forest = T)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 21
Mean of squared residuals: 173.1405
% Var explained: 88.71
train.set %>%
filter(Compartment == "Endosphere") %>%
.$rf
[[1]]
Call:
randomForest(x = imp.predictors, y = response, importance = F, keep.forest = T)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 21
Mean of squared residuals: 149.3839
% Var explained: 90.26
Use the whole dataset to get the predicted age of the plants
get_predict <- function(model, data) {
predict(model, data)
}
test.set <- whole.set %>%
inner_join(select(train.set, Compartment, rf), by = "Compartment") %>%
mutate(prediction = map2(rf, otu, get_predict))
rs.pred <- test.set %>% unnest(prediction, age) %>% filter(Treatment == "WC") %>% filter(Compartment == "Rhizosphere") %>% .$prediction
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(prediction, age))`, with `mutate()` if needed
rs.age <- test.set %>% unnest(prediction, age) %>% filter(Treatment == "WC") %>% filter(Compartment == "Rhizosphere") %>% .$age
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(prediction, age))`, with `mutate()` if needed
cor(rs.pred, rs.age)
[1] 0.9638256
es.pred <- test.set %>% unnest(prediction, age) %>% filter(Treatment == "WC") %>% filter(Compartment == "Endosphere") %>% .$prediction
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(prediction, age))`, with `mutate()` if needed
es.age <- test.set %>% unnest(prediction, age) %>% filter(Treatment == "WC") %>% filter(Compartment == "Endosphere") %>% .$age
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(prediction, age))`, with `mutate()` if needed
cor(es.pred, es.age)
[1] 0.9774801
test.set %>%
unnest(prediction, age) %>%
ggplot(aes(age, prediction)) +
geom_point(aes(shape = Set, color = Treatment), alpha = 0.7) +
geom_abline(linetype = 2) +
geom_smooth(aes(color = Treatment), linetype = 2, stat = "smooth", se = F) +
scale_color_manual(values = trt.pal) +
facet_grid(Treatment + Set ~ Compartment) +
theme_light() +
theme(text = element_text(size = 20))
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(prediction, age))`, with `mutate()` if needed
predictions <- test.set %>%
unnest(prediction, age)
unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(prediction, age))`, with `mutate()` if needed
predictions.wc <- filter(predictions, Treatment == "WC" & Set == "Test") %>% select(-Treatment) %>%
ungroup()
Adding missing grouping variables: `Treatment`
predictions.wc.plotting <- rbind(select(predictions.wc, -Treatment) %>% mutate(Treatment = "DS1"),
select(predictions.wc, -Treatment) %>% mutate(Treatment = "DS2"),
select(predictions.wc, -Treatment) %>% mutate(Treatment = "DS3"))
predictions %>%
filter(Treatment != "WC") %>%
mutate(Treatment = str_replace(Treatment, "D", "DS")) %>%
mutate(Treatment = fct_relevel(Treatment, "WC")) %>%
ggplot(aes(age, prediction)) +
geom_vline(data = trt.lines, aes(xintercept = IniTreatment2), linetype = 3) +
geom_vline(data = trt.lines, aes(xintercept = EndTreatment2), linetype = 3) +
geom_point(aes(color = Treatment), size = 2, shape = 1, alpha = 1) +
geom_smooth(data = predictions.wc.plotting, linetype = 2, size = 1, stat = "smooth", se = F, color = trt.pal[1]) +
xlab("Plant Age (days)") +
ylab("Microbiome Age\n(days)") +
scale_color_manual(values = trt.pal[-1]) +
facet_grid(Compartment ~ Treatment) +
theme_classic() +
theme(text = element_text(size = 20),
legend.position = "none")
Unknown levels in `f`: WCUnknown levels in `f`: WCUnknown levels in `f`: WCUnknown levels in `f`: WCUnknown levels in `f`: WCUnknown levels in `f`: WC
loess.pred <- rbind(
data.frame(loess = predict(loess(prediction~age,filter(predictions.wc, Compartment == "Rhizosphere")), unique(predictions.wc$age)),
age = unique(predictions.wc$age),
Compartment = "Rhizosphere"),
data.frame(loess = predict(loess(prediction~age,filter(predictions.wc, Compartment == "Endosphere")), unique(predictions.wc$age)),
age = unique(predictions.wc$age),
Compartment = "Endosphere"))
predictions.wc %>% ggplot(aes(age, prediction, color = Compartment)) +
geom_smooth(linetype = 2, size = 1, stat = "smooth", se = F) +
geom_point(data = loess.pred, aes(age, loess))
NA
NA
library(broom)
library(contrast)
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
time.age <- map %>%
group_by(Time, Age) %>%
count() %>%
select(-n)
predictions <- inner_join(predictions, time.age, by = c("age" = "Age"))
predictions$TimeFctr <- as.factor(predictions$Time)
run_lm <- function(df) {
lm(prediction ~ Treatment * TimeFctr, data = df)
}
run_anova <- function(fit) {
anova(fit) %>% tidy()
}
get_contrasts <- function(fit) {
contrasts <- vector(mode = "list")
for(i in 1:13) {
i <- as.character(i)
print(i)
contrasts[[paste("D1", i, sep = ".")]] <- contrast(fit, list(Treatment = "WC", TimeFctr = i), list(Treatment = "D1", TimeFctr = i))[1:8] %>% as.data.frame()
contrasts[[paste("D2", i, sep = ".")]] <- contrast(fit, list(Treatment = "WC", TimeFctr = i), list(Treatment = "D2", TimeFctr = i))[1:8] %>% as.data.frame()
contrasts[[paste("D3", i, sep = ".")]] <- contrast(fit, list(Treatment = "WC", TimeFctr = i), list(Treatment = "D3", TimeFctr = i))[1:8] %>% as.data.frame()
}
contrasts <- plyr::ldply(contrasts, function(x) x)
contrasts <- contrasts %>%
separate(.id, c("Treatment", "Time"))
}
predictions.nested <-predictions %>%
group_by(Compartment) %>%
nest() %>%
mutate(lm = map(data, run_lm),
anova = map(lm, run_anova),
contrasts = map(lm, get_contrasts))
[1] "1"
[1] "2"
[1] "3"
[1] "4"
[1] "5"
[1] "6"
[1] "7"
[1] "8"
[1] "9"
[1] "10"
[1] "11"
[1] "12"
[1] "13"
[1] "1"
[1] "2"
[1] "3"
[1] "4"
[1] "5"
[1] "6"
[1] "7"
[1] "8"
[1] "9"
[1] "10"
[1] "11"
[1] "12"
[1] "13"
rf.sig <- predictions.nested %>%
unnest(contrasts) %>%
group_by(Compartment, Treatment) %>%
mutate(p.adj = p.adjust(Pvalue, method = "fdr")) %>%
mutate(Time = as.integer(Time)) %>%
inner_join(time.age, by = "Time")
a <- predictions %>%
filter(Treatment != "WC") %>%
mutate(Compartment = fct_relevel(Compartment, "RS", "ES")) %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
mutate(Treatment = str_replace(Treatment, "D", "DS")) %>%
ggplot(aes(age, prediction)) +
geom_vline(data = trt.lines, aes(xintercept = IniTreatment2), linetype = 3) +
geom_vline(data = trt.lines, aes(xintercept = EndTreatment2), linetype = 3) +
geom_point(aes(color = Treatment), size = 2, shape = 1, alpha = 1) +
geom_smooth(data = predictions.wc.plotting, linetype = 2, size = 1, stat = "smooth", se = F, color = trt.pal[1]) +
xlab("Plant Age (days)") +
ylab("Microbiome Age\n(days)") +
scale_color_manual(values = trt.pal[-1]) +
facet_grid(Compartment ~ Treatment) +
theme_classic() +
theme(text = element_text(size = 15),
legend.position = "none")
Unknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: RS, ES
a
b <- predictions %>%
filter(Treatment != "WC") %>%
inner_join(loess.pred, by = c("Compartment", "age")) %>%
mutate(RelMat = (prediction - loess)) %>%
group_by(Compartment, Treatment, age) %>%
summarise(MeanMat = mean(RelMat),
SDMat = sd(RelMat)) %>%
mutate(SEMat = SDMat/2) %>%
inner_join(rf.sig, by = c("Compartment","Treatment", "age" = "Age")) %>%
ungroup() %>%
mutate(FDR = ifelse(p.adj < 0.05, "< 0.05", "≥ 0.05")) %>%
mutate(Compartment = fct_relevel(Compartment, "RS", "ES")) %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
mutate(Treatment = str_replace(Treatment, "D", "DS")) %>%
mutate(Treatment = fct_relevel(Treatment, "WC")) %>%
ggplot(aes(age, MeanMat, color = Treatment, shape = FDR)) +
geom_point(size = 2) +
geom_errorbar(aes(ymin = MeanMat - SEMat, ymax = MeanMat + SEMat)) +
geom_vline(data = trt.lines, aes(xintercept = IniTreatment2), linetype = 3) +
geom_vline(data = trt.lines, aes(xintercept = EndTreatment2), linetype = 3) +
ylab("Relative Microbiome\nMaturity (days)") +
xlab("Plant Age (days)") +
facet_grid(Compartment ~ Treatment) +
scale_color_manual(values = trt.pal[-1]) +
scale_shape_manual(values = rev(ws.shp)) +
theme_classic() +
theme(text = element_text(size = 15),
legend.position = "bottom") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
shape = guide_legend(title.position = "top", title.hjust = 0.5))
`summarise()` regrouping output by 'Compartment', 'Treatment' (override with `.groups` argument)
Unknown levels in `f`: RS, ESUnknown levels in `f`: RS, ESUnknown levels in `f`: WC
b
c <- otu.tidy %>%
mutate(Count = Count/100) %>%
inner_join(map, by = "SampleID") %>%
filter(Treatment2 != "WC_TRN") %>%
inner_join(all.ord, by = c("Compartment", "OTU_ID")) %>%
group_by(Compartment, Treatment, Age, Time, ZS_cluster, Trend, SampleID) %>%
summarise(Total = sum(Count)) %>%
group_by(Compartment, Treatment, Age, Time, ZS_cluster, Trend) %>%
mutate(MeanAb = mean(Total)) %>%
ungroup() %>%
mutate(Trend = fct_relevel(Trend, "Complex", after = Inf)) %>%
mutate(Compartment = fct_relevel(Compartment, "RS", "ES")) %>%
mutate(Compartment = fct_recode(Compartment, Rhizosphere = "RS", Endosphere = "ES")) %>%
mutate(Treatment = str_replace(Treatment, "D", "DS")) %>%
mutate(Treatment = fct_relevel(Treatment, "WC")) %>%
ggplot(aes(Age, Total, fill = Trend, color = Trend)) +
geom_point(shape = 1, size = 2) +
geom_smooth(se = F) +
geom_vline(data = trt.lines, aes(xintercept = IniTreatment2), linetype = 3) +
geom_vline(data = trt.lines, aes(xintercept = EndTreatment2), linetype = 3) +
ylab("Relative Abundance") +
xlab("Plant Age (days)") +
scale_fill_manual(values = c("dodgerblue", "gold", "grey50")) +
scale_color_manual(values = c("dodgerblue", "gold", "grey50")) +
facet_grid(Treatment ~ Compartment) +
theme_classic() +
theme(text = element_text(size = 15),
legend.position = "bottom") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
fill = guide_legend(title.position = "top", title.hjust = 0.5))
`summarise()` regrouping output by 'Compartment', 'Treatment', 'Age', 'Time', 'ZS_cluster', 'Trend' (override with `.groups` argument)
c
left <- plot_grid(a,b + theme(legend.position = "none"), get_legend(b),
nrow = 3,
rel_heights = c(4,4,1),
labels = c("A", "B"),
label_size = 20)
conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <e2>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <89>conversion failure on '≥ 0.05' in 'mbcsToSbcs': dot substituted for <a5>`geom_smooth()` using method = 'loess' and formula 'y ~ x'
right <- plot_grid(c + theme(legend.position = "none"), get_legend(c),
nrow = 2,
rel_heights = c(8,1),
labels = "C",
label_size = 20)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
##1011:849
plot_grid(left, right,
rel_widths = c(3,2))